function [theta , opts] = ML_optimize_global(theta, data, opts)
% theta = ML_optimize_global(theta, data, opts)
% Optimizes all of theta in one global optimization.
% The assumed distribution is specified by opts.logPDFfun
%
% Created by Dirk Poot, Erasmus MC, 22-3-2011
opt = optimset('fminunc');
opt = optimset(opt,'gradObj','on','Hessian','on', 'Display','iter', 'MaxIter',opts.maxIter,'TolFun',opts.tolFun,'largescale','on','PrecondBandWidth',inf,'OutputFcn',@progbupd);

opts.explicitHessian = false; 
        % Never compute hessian explicitly, it might be very large (and if it's small, thats probably for testing).
if opts.doRegularize
    regularizer.fun = opts.spatialRegularizer{1};
    regularizer.explicitHessian = opts.explicitHessian;
    regularizer.hessMulFun = opts.spatialRegularizer{2};
end;

if ~opts.explicitHessian
    % if not explicit hessian, I use the following convention for hessinfo:
    %
    % hessinfo  = the third output of a function to be optimized
    %
    % [Hy] = hessinfo.hessMulFun( hessinfo, Y )  multiplies the (possibly multi column vector) Y with the hessian
    %  => Hy = H * Y
    % [R, fun] = hessinfo.makePreconditioner( hessinfo , upperbandw, DM, DG )
    % prepare R as preconditioner of 
    %    M = DM*H*DM + DG
    % so 
    %    fun(x, R)   =approx=   inv(M) * x
    maxPCGiter = min(prod(opts.blockSize)*size(theta,1),60+20*ndims(theta));
    hessmulfun = @(hessinfo, Y) fullThetaCostfun_HessMul(hessinfo, Y,opts); 
    opt = optimset(opt,'HessMult',hessmulfun, 'MaxPCGIter',maxPCGiter,'PrecondBandWidth',0); 
        % set precondbandwidth to less than inf, to display cgiter correctly.
else
    hessmulfun = [];
end;
if ~isempty(opts.constraints)
    %opts.constraints_A,  opts.constraints_b
    npar = size(theta,1);
    constrFun = @(x) constraintBlockFun(x, opts.constraints , npar );
end;
% select 'fields' (i.e. additional parameters per voxel)
if ~isempty(opts.fields)
    opts.function = make1arg_anonfun( opts.function, opts.fields(:,:) );
end;

optfun = @(t)  fullThetaCostfun(t, opts);
progressbar('start',[0 maxPCGiter;0 opts.maxIter]);
if isempty(opts.constraints)
    [theta_opt, fval, exflag, outp, grad, hess] = fminunc(optfun, theta, opt);
else
    [theta_opt, fval, exflag, outp, lambd, grad, hess] = fmincon(optfun, theta, [],[],[],[],[],[], constrFun ,opt);
end;
theta = theta_opt;
progressbar('ready');
opts.curpos.LLdata = hess.LLfun;
opts.curpos.LLregularization = hess.LLregularization;
if opts.doComputeDerivativeRegularizationScale
    % Compute dThetadLamba = - H^(-1) * dRegularizationdTheta
    % 
    opts.curpos.dRegularizationdTheta = hess.dRegularizationdTheta;
%     opts.curpos.dThetadLamba = pcg( @(Y) fullThetaCostfun_HessMul(hess, Y,opts), hess.dRegularizationdTheta, 1e-4, 100);
    opts.curpos.dThetadLamba =  - .5 * cgiterLS( @(x) 0, @(x) x, hess.dRegularizationdTheta,[],@(Y) fullThetaCostfun_HessMul(hess, Y,opts), 100);
end;

function [stop] =progbupd(xOutputfcn,optimValues,state)
progressbar( [optimValues.cgiterations ; optimValues.funccount ] );
stop = false;


